/*
 *	$Id: matematik.h 157 2009-06-19 15:05:08Z wabsie $
 */

#ifndef _MATEMATIK_H
#define	_MATEMATIK_H

#include <cmath>		// STL udgave af math.h
#include <string>		// STL string
#include <vector>		// STL vector

using namespace std;	// Brug STL namespace

typedef double skalar;	// Gør det nemt at ændre beregningspræcsionen

struct punkt;
class plotter;			// Prototyper på klasserne
template <int N, int M, bool D> class base;
template <int N, int M>			class  matrix;
template <int N>				class  base23;
template <int N>				class   pixel;
template <int N>				class   coord;
template <int N>				class   vektor;

#ifdef M_PI
static const skalar pi = M_PI;
#else
static const skalar pi = 4*atan(1); // Beregn pi nøjagtigt 1. gang programmet starter
#endif

static const skalar INF = numeric_limits<skalar>::infinity();	// +Uendelig
static const skalar EPS = numeric_limits<skalar>::epsilon();	// Mindste maskinrepræsenterbare tal

#define AFRUND(x)		(int)round(x)
#define	POTENS10(x)		((x) == 0) ? 0 : floor(log10(abs((x))))


template <int N, int M, bool dynamic = false>
struct trait
{
	typedef skalar foo[N][M];

	static void alloc(foo& data)
	{
		// Statisk allokerede arrays skal ikke allokeres dynamisk
	}

	static void zero(foo& data)
	{
		for (int i = 0; i < N; i++)
		for (int j = 0; j < M; j++)
			data[i][j] = 0;
	}
};

template <int N, int M>
struct trait<N, M, true>
{
	typedef vector< vector<skalar> > foo;

	static void alloc(foo& data, int nn = N, int mm = M)
	{
		data.resize(nn);
		for (int i = 0; i < nn; i++)
			data[i].resize(mm);
	}

	static void zero(foo& data, int nn = N, int mm = M, skalar initval = 0)
	{
		for (int i = 0; i < nn; i++)
		for (int j = 0; j < mm; j++)
			data[i][j] = initval;
	}
};


template <int N, int M = 1, bool D = false>
class base
{	/*
	 * Default til statisk allokeret søjlevektor.
	 */

private:
	/*
	 * Jvf. http://womble.decadentplace.org.uk/c++/template-faq.html#type-syntax-error
	 * om hvorfor 'typename' indgår i vores typedefine her under.
	 */
	typedef typename trait<N,M,D>::foo dataType;

public:
	int n, m;
	dataType data;

	base() :
		n(N), m(M)
	{	/* Nulstil her; umuligt i init-list med arrays */
		trait<N,M,D>::alloc(this->data);	// Allokér om nødvendigt
		trait<N,M,D>::zero(this->data);		// Nulstil
	}

	base(const base<N,M>& arg) :
		n(N), m(M)
	{	/* Copy constructor */
		trait<N,M,D>::alloc(this->data);	// Allokér hvis nødvendigt

		for (int i = 0; i < n; i++)			// Deep copy
		for (int j = 0; j < m; j++)
			this->data[i][j] = arg.data[i][j];
	}

	const string toString(const char* fmt = "%fe%02.0f", int j = 0) const
	{
		string res = "";
		char buf[BUFSIZ] = "";

		for (int i = 0; i < n; i++) {
			snprintf(buf, sizeof(buf), "%s%s",	fmtScientific(data[i][j], fmt).c_str(),
												(i == n-1) ? "" : ";");
			res += buf;
		}

		return res;
	}

	inline base<N,M> operator+(const base<N,M>& rhs) const
	{	/* Binær plus */
		base<N,M> res = *this;

		for (int i = 0; i < n; i++)
		for (int j = 0; j < m; j++)
			res.data[i][j] += rhs.data[i][j];

		return res;
	}

	inline base<N,M> operator-(const base<N,M>& rhs) const
	{	/* Binær minus */
		base<N,M> res = *this;

		for (int i = 0; i < n; i++)
		for (int j = 0; j < m; j++)
			res.data[i][j] -= rhs.data[i][j];

		return res;
	}

	inline base<N,M> operator-(void) const
	{	/* Unær minus */
		return this->operator *(-1);
	}

	inline base<N,M> operator+(void) const
	{	/* Unær plus */
		return *this; // Ingen ændring fortaget
	}

	inline base<N,M> operator*(const skalar rhs) const
	{	/* Gange med skalar */
		base<N,M> res = *this;

		for (int i = 0; i < n; i++)
		for (int j = 0; j < m; j++)
			res.data[i][j] *= rhs;

		return res;
	}

	inline base<N,M> operator/(const skalar rhs) const
	{	/* Gange med skalar */
		base<N,M> res = *this;

		for (int i = 0; i < n; i++)
		for (int j = 0; j < m; j++)
			res.data[i][j] /= rhs;

		return res;
	}

	inline skalar norm(int j = 0) const
	{	/* Finder 2-normen af vores objekt (betragtet som vektor) */
		skalar res = 0;

		for (int i = 0; i < n; i++)
			res += this->data[i][j] * this->data[i][j]; // Undgå kald til 'pow'

		return sqrt(res);
	}

	inline base<N,M> mul(const base<N,M>& rhs) const
	{	/* Gang elementvist */
		base<N,M> res = *this;

		for (int i = 0; i < n; i++)
		for (int j = 0; j < m; j++)
			res.data[i][j] *= rhs.data[i][j];

		return res;
	}

	inline base<N,M> div(const base<N,M>& rhs) const
	{	/* Dividér elementvist */
		base<N,M> res = *this;

		for (int i = 0; i < n; i++)
		for (int j = 0; j < m; j++)
			res.data[i][j] /= rhs.data[i][j];

		return res;
	}

	inline base<N,M>& assign(const base<N,M>& rhs)
	{
		if (this == &rhs)
			return *this;

		for (int i = 0; i < n; i++)
		for (int j = 0; j < m; j++)
			this->data[i][j] = rhs.data[i][j];

		return *this;
	}
};


template <int N = 0, int M = 0>
class matrix : public base<N,M>
{
public:
	matrix() : base<N,M>() {}							// Zero-initialisér
	matrix(const base<N,M>& arg) : base<N,M>(arg) {}	// Kopiér indhold
};

template <>
class matrix<0,0> : public base<1,1,true>
{
public:
	matrix(int _n, int _m, skalar initval = 0) :
		base<1,1,true>()
	{
		trait<1,1,true>::alloc(this->data, _n, _m);	// Allokér; vi er dynamiske nu
		trait<1,1,true>::zero(this->data, _n, _m, initval);	// Nulstil

		this->n = _n;
		this->m = _m;
	}
};


template <int N>
class base23 : public base<N>
{
	base23() : base<N>() {}							// Zero-initialisér
	base23(const base<N>& arg) : base<N>(arg) {}	// Kopiér indhold
};

template <>
class base23<2> : public base<2>
{
public:
	skalar& x; skalar& y;

	base23() :
		base<2>(),		// Zero-initialisér
		x(this->data[0][0]), y(this->data[1][0]) // Sæt alias' op
	{
		// Intet her
	}

	// Vores egen copy constructor
	base23(const base<2>& arg) :
		base<2>(arg),	// Kopiér indhold
		x(this->data[0][0]), y(this->data[1][0]) // Sæt alias' op
	{
		// Intet her
	}

	base23<2>& operator=(const base23<2>& rhs)
	{	/* Assignment operator */
		this->assign(rhs);
		return *this;
	}
};

template <>
class base23<3> : public base<3>
{
public:
	/*
	 * Undgå at skrive this->data[][]: Tillader blot data[][].
	 * Se også http://womble.decadentplace.org.uk/c++/template-faq.html#base-lookup
	 */
	using base<3>::data;

	skalar& x; skalar& y; skalar& z;

	base23() :
		base<3>(),		// Zero-initialisér
		x(data[0][0]), y(data[1][0]), z(data[2][0]) // Setup alias' (referancer)
	{
		// Intet her
	}

	// Vores egen copy constructor
	base23(const base<3>& arg) :
		base<3>(arg),	// Kopiér indhold
		x(data[0][0]), y(data[1][0]), z(data[2][0]) // Setup alias'
	{
		// Intet her
	}

	base23<3>& operator=(const base23<3> rhs)
	{	/* Assignment operator */
		this->assign(rhs);
		return *this;
	}
};


template <int N>
class coord : public base23<N>
{
	coord() : base23<N>() {}						// Default constructor: Origo
	coord(const base<N>& arg) : base23<N>(arg) {}	// Copy constructor: Betragt arg som et koordinatpunkt
	coord(const coord<N>& arg) : base23<N>( (base<N>)arg ) {}
};

template <>
class coord<2> : public base23<2>
{
public:
	coord() : base23<2>() {}
	coord(const base<2>& arg) : base23<2>(arg) {}
	coord(const coord<2>& arg) : base23<2>( (base<2>)arg ) {}
	coord(const skalar _x, const skalar _y)
	{
		this->x = _x; // Ækvivalent med: this->data[0][0] = _x;
		this->y = _y; // Ækvivalent med: this->data[1][0] = _y;
	}

	inline pixel<2> toPixel(plotter *p) const;	// jvf. plotter.h
};

template <>
class coord<3> : public base23<3>
{
public:
	coord() : base23<3>() {}
	coord(const base<3>& arg) : base23<3>(arg) {}
	coord(const coord<3>& arg) : base23<3>( (base<3>)arg ) {}
	coord(const skalar _x, const skalar _y, const skalar _z)
	{
		this->x = _x; // Ækvivalent med: this->data[0][0] = _x;
		this->y = _y; // Ækvivalent med: this->data[1][0] = _y;
		this->z = _z; // Ækvivalent med: this->data[2][0] = _z;
	}

	inline punkt toPunkt(plotter *p) const;	// jvf. plotter.h

	inline coord<3> rotx(skalar rad) const
	{
		coord<3> Pm;
		skalar cosv = cos(rad);
		skalar sinv = sin(rad);

		Pm.x = this->x;
		Pm.y = this->y * cosv - this->z * sinv;
		Pm.z = this->z * cosv + this->y * sinv;

		return Pm;
	}

	inline coord<3> roty(skalar rad) const
	{
		coord<3> Pm;
		skalar cosv = cos(rad);
		skalar sinv = sin(rad);

		Pm.x = this->x * cosv + this->z * sinv;
		Pm.y = this->y;
		Pm.z = this->z * cosv - this->x * sinv;

		return Pm;
	}

	inline coord<3> rotz(skalar rad) const
	{
		coord<3> Pm;
		skalar cosv = cos(rad);
		skalar sinv = sin(rad);

		Pm.x = this->x * cosv - this->y * sinv;
		Pm.y = this->y * cosv + this->x * sinv;
		Pm.z = this->z;

		return Pm;
	}
};


template <int N>
class pixel : public base23<N>
{
	pixel() : base23<N>() {}						// Default constructor: "Origo"-pixlen
	pixel(const base<N>& arg) : base23<N>(arg) {}	// Copy constructor: Betragt arg som en pixel
	pixel(const pixel<N>& arg) : base23<N>( (base<N>)arg ) {}
};

template <>
class pixel<2> : public base23<2>
{
public:
	pixel() : base23<2>() {}
	pixel(const base<2>& arg) : base23<2>(arg) {}
	pixel(const pixel<2>& arg) : base23<2>( (base<2>)arg ) {}
	pixel(const skalar _x, const skalar _y)
	{
		this->x = _x; // Ækvivalent med: this->data[0][0] = _x;
		this->y = _y; // Ækvivalent med: this->data[1][0] = _y;
	}

	inline coord<2> toCoord(plotter *p) const;	// jvf. plotter.h
};

template <>
class pixel<3> : public base23<3>
{
public:
	pixel() : base23<3>() {}
	pixel(const base<3>& arg) : base23<3>(arg) {}
	pixel(const pixel<3>& arg) : base23<3>( (base<3>)arg ) {}
	pixel(const skalar _x, const skalar _y, const skalar _z)
	{
		this->x = _x; // Ækvivalent med: this->data[0][0] = _x;
		this->y = _y; // Ækvivalent med: this->data[1][0] = _y;
		this->z = _z; // Ækvivalent med: this->data[2][0] = _z;
	}
};


template <int N>
class vektor : public base23<N>
{
public:
	vektor() : base23<N>() {}						// Default constructor: Nulvektoren
	vektor(const base<N>& arg) : base23<N>(arg) {}	// Copy constructor: Betragt arg som vektor
	vektor(const vektor<N>& arg) : base23<N>( (base<N>)arg ) {} // Undgå referancekopiering og def ctor
	vektor(const base<N>& p1, const base<N>& p2)	// Vektoren mellem to punkter (ikke stedvektor)
	{	/* Opret vektoren der går fra p1 til p2 */
		*this = p2 - p1;
	}
};

template <>
class vektor<2> : public base23<2>
{
public:
	vektor() : base23<2>() {}
	vektor(const base<2>& arg) : base23<2>(arg) {}
	vektor(const vektor<2>& arg) : base23<2>( (base<2>)arg ) {} // Undgå referancekopiering
	vektor(const skalar _x, const skalar _y)
	{	/* Opret vektoren <_x, _y> */
		this->x = _x; // Ækvivalent med: this->data[0][0] = _x;
		this->y = _y; // Ækvivalent med: this->data[1][0] = _y;
	}
	vektor(const base<2>& p1, const base<2>& p2)
	{	/* Opret vektoren der går fra p1 til p2 */
		*this = p2 - p1;
	}

	vektor<2> orto() const 
	{	/* Finder tværvektoren til this */
		return vektor<2>(-this->y, this->x);
	}
};

template <>
class vektor<3> : public base23<3>
{
public:
	vektor() : base23<3>() {}
	vektor(const base<3>& arg) : base23<3>(arg) {}
	vektor(const vektor<3>& arg) : base23<3>( (base<3>)arg ) {} // Undgå referancekopiering
	vektor(const skalar _x, const skalar _y, const skalar _z)
	{	/* Opret vektoren <_x, _y, _z> */
		this->x = _x; // Ækvivalent med: this->data[0][0] = _x;
		this->y = _y; // Ækvivalent med: this->data[1][0] = _y;
		this->z = _z; // Ækvivalent med: this->data[2][0] = _z;
	}
	vektor(const base<3>& p1, const base<3>& p2)
	{	/* Opret vektoren der går fra p1 til p2 */
		*this = p2 - p1;
	}

	vektor<3> kryds(const vektor<3>& rhs) const
	{	/* Krydsproduktet */
		vektor <3> tmp;
		tmp.x = this->y * rhs.z - this->z * rhs.y; // (a2*b3 - a3*b2)
		tmp.y = this->z * rhs.x - this->x * rhs.z; // (a3*b1 - a1*b3)
		tmp.z = this->x * rhs.y - this->y * rhs.x; // (a1*b2 - a2*b1)

		return tmp;
	}

	skalar dot(const vektor<3>& rhs) const
	{	/* Prikproduktet */
		return this->x*rhs.x + this->y*rhs.y + this->z*rhs.z;
	}
};

struct fieldPoint
{
	enum felttype {dont_know, skalar_felt, vektor_felt} type;
	skalar sfelt;
	vektor<3> vfelt;

	fieldPoint() :
		// Initialisér feltet til Intet
		sfelt(0), vfelt(vektor<3>()), type(dont_know)
	{
		// Intet her
	}

	fieldPoint(felttype _type) :
		sfelt(0), vfelt(vektor<3>()), type(_type)
	{
		// Intet her
	}

	inline fieldPoint operator+(const fieldPoint& rhs) const
	{	/* Binær plus: Læg to felter af samme type sammen */
		fieldPoint tmp;

		tmp.type = (rhs.type == dont_know) ? this->type : rhs.type;
		assert(tmp.type != dont_know); // Vær sikker på at vi kender feltets type nu

		if (tmp.type == skalar_felt)
			tmp.sfelt = this->sfelt + rhs.sfelt;
		else
			tmp.vfelt = this->vfelt + rhs.vfelt;

		return tmp;
	}

	inline skalar getItensity() const
	{
		assert(this->type != dont_know); // Vær sikker på at vi kender feltets type

		if (this->type == skalar_felt)
			return abs(this->sfelt);	// abs: Kan være negativ, og vi vil kun have intensiteten
		else
			return this->vfelt.norm();	// Aldrig negativ; abs kan undværes
	}
};


static const string fmtScientific(skalar tal, const char* fmt = "%fe%02.0f")
{
	char buf[BUFSIZ] = "";
	skalar potens = POTENS10(tal);
	snprintf(buf, sizeof(buf), fmt, tal/pow(10, potens), potens);

	return (string)buf;
}


struct punkt {
	string id;
	coord<3> orig;		// Det rigtige 3D punkt.
	coord<3> offrot;	// 'orig' translateret og roteret, men ikke projiceret
	coord<3> proj;		// Projektionen af 'offrot' (2D projektionen er blot [proj.x, proj.y])
	pixel<2> pixl;		// Tilsvarende pixel-koordinat for [proj.x, proj.y]
};


#endif	/* _MATEMATIK_H */
